Spatial and Temporal Heterogeneity of 
Host-Parasitoid Interactions in Lupine 

Habitat 

By 

Roy Werner Wright 
B.S. (University of California, Irvine) 2004 

00 
O 

^ DISSERTATION 

^ Submitted in partial satisfaction of the requirements for the degree of 

< 

DOCTOR OF PHILOSOPHY 
in 

Ph applied mathematics 

d 

•pl in the 

S OFFICE OF GRADUATE STUDIES 

^ of the 

00 

00 UNIVERSITY OF CALIFORNIA 

a^ 

^ DAVIS 
oo 

o 

00 

O Approved: 
> 



Alan Hastings (Chair) 



John Hunter 



Alex Mogilner 
Committee in Charge 
2008 



i 



Contents 



1 Introduction: The Lupine Habitat [T] 

2 Spontaneous Patchiness in an Integrodifference Model |4] 

2.1 Introduction H] 

2.2 The Model E 

2.3 Singular Perturbation Analysis of Steady State [7] 

2.3.1 Regular Perturbation [7] 

2.3.2 Transition Layers [H] 

2.3.3 Summary of Results [10] 

2.4 Numerical Experiments [12] 

2.4.1 Laplace Kernel [H 

2.4.2 Other Kernels M 

2.5 Routes to Heterogeneity [13 

2.5.1 The Nonspatial Model [I5] 

2.5.2 Evolution to the Steady State [16] 

2.5.3 Dependence on Parameters and Initial Conditions [T7] 

2.5.4 Linear Stability M 

2.5.5 A Bound on Patch Radius [23] 

ii 



2.6 Discussion [25] 

3 The Effects of Habitat Quality on Patch Formation [29] 

3.1 Introduction [29] 

3.2 An Allee Effect Induced by Saturating Predation [HO] 

3.2.1 Modeling the Growth of the Host M 

3.2.2 Escape Probability [32] 

3.2.3 Continuous Predation [33] 

3.2.4 Properties of the Maps [3l 

3.3 Consequences of Habitat Quality [3H] 

3.3.1 Introduction [38] 

3.3.2 Dependence on Carrying Capacity [39] 

3.3.3 Numerical Examples [10] 

3.3.4 The Paradox of Enrichment [13] 

3.4 Consequences of Habitat Heterogeneities [S] 

3.4.1 Introduction [H 

3.4.2 Patch Formation [H 

3.4.3 Numerical Examples [IH] 

3.5 Extension to Two Dimensions [19] 

3.6 Discussion 

4 Persistence of Nematodes Augmented by an Alternate Host [58] 

4.1 Introduction [5H] 

4.2 Details from the Original Model [59] 

4.2.1 In- Year Dynamics [60] 

4.2.2 Bifurcation Diagram of the Map [HI] 

iii 



4.3 The New Model ED 

4.3.1 Unsuitability of Deterministic Modeling [HI] 

4.3.2 The Model and Its Simulation [63] 

4.4 Means of Persistence [M] 

4.4.1 Mitigated Crashing El 

4.4.2 Transient Survival [65] 

4.4.3 Deterministic Stability [HZ] 

4.5 Discussion [69] 

References [7T] 



iv 



Roy Werner Wright 
September 2008 
Applied Mathematics 

Spatial and Temporal Heterogeneity of Host-Parasitoid Interactions 

in Lupine Habitat 
Abstract 

The inhabitants of the bush lupine in coastal California have been the subject 
of scientific scrutiny in recent years. Observations of a host-par asitoid interaction in 
the shrub's foliage, in which victims are significantly less motile than their exploiters, 
record stable spatial patterns in a fairly homogeneous environment. Though such 
pattern formation has been found in reaction-diffusion models, the correspondence of 
these models to continuous-time predator-prey interactions does not reflect the reality 
of the system being studied. Near the root of the lupine, another host-parasitoid 
interaction is also of considerable interest. In some cases this interaction, which 
promotes the health of the lupine, has been observed to be much more persistent 
than suggested by mathematical models. 

In this work a discrete-time spatial model of the first host-parasitoid system is 
introduced. We analyze the model, describing its transient behavior and finding the 
conditions under which spatial patterns occur, as well as an estimate of outbreak 
size under those conditions. We consider the feasibility of the necessary conditions 
in the natural system by modeling the mechanisms responsible for them, and discuss 
the effects of variable habitat on pattern formation. We also explore one possible 
explanation for the persistence of the second host-parasitoid system - the existence of 
an alternate host. Under certain surprising conditions and by means quite different 
from previous models of similar situations, an alternate host can greatly enhance 
persistence of the nematode par asitoid. 
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Chapter 1 

Introduction: The Lupine Habitat 

The ecological system motivating the present work centers on bush lupines on the 
coast of California. Above ground, these shrubs are home to flightless tussock moths, 
which feed on the bush but have little or no effect on its year-to-year dynamics. In 
turn, the tussock moths are parasitized by a variety of much more motile enemies [19] . 
Below ground, the bush lupine's roots provide shelter and sustenance to destructive 
ghost moth larvae, which are in turn parasitized by entomopathogenic nematodes. In 
this way, the nematodes promote the health of the shrubs through a trophic cascade. 

Let us briefly review some pertinent results, to which we will return in the follow- 
ing chapters. In |20j it is shown experimentally that the pupae of the tussock moth 
are subject to mortality by a predatory ant with a saturating functional response. 
In [22], measurements of tussock moth density from the field are compared with re- 
sults from a reaction-diffusion model of a victim and its more motile exploiter. The 
model gives a good qualitative prediction of the somewhat counterintuitive spatial 
distribution of the tussock moth. An integrodifference model for the tussock moth is 
analyzed in |18] , but it proves to be incapable of some of the spatial patterns seen in 
the field. 
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In [9] and [10], deterministic and stochastic models, respectively, show the ten- 
dency of soil-dwelling nematode populations to oscillate wildly from year to year 
between very high and dangerously low levels. Some experimental work in [TU] sug- 
gests that model values for nematode mortality have been too severe in some cases, 
but it is also noted that many other possible explanations for the observed persistence 
of nematodes in some rhizospheres have been all but ruled out. 

In [19], the lupine-centered community is discussed in detail and a preliminary 
step is taken to link the above- and below-ground subsystems. This step consists of 
analyzing the effects of changing carrying capacity - which represents the health of 
the bush lupine - on the tussock moth. 




Figure 1.0.1: Species of the lupine habitat: soil-dwelling Nematodes with their larval 
Ghost moth hosts and possible alternate Hosts, Tussock moths with their motile 
Parasitoids and predatory Ants, and the bush Lupine. The meaning of arrows here 
is nonstandard (see text); trophic level is represented by vertical position. 



A graphical overview of the community described above, as modeled in this dis- 
sertation, is given in Figure [T.0.1 Arrows point in the direction of modeled ecological 
influences. For example, lupine (L) abundance affects the dynamics of the tussock 
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moth (T), but defohation by the tussock moth has no effect on the health of the 
lupine. The tussock moth and its parasitoids (P) are mutually interacting, while 
it is assumed that the generalist ant predator (A) is not limited by tussock moth 
abundance. We analyze an integrodifference model for the interaction of the tussock 
moth with its parasitoid enemies in Chapter |2} In Chapter [3} we consider in greater 
detail the influence of predatory ants and the health of the lupine habitat. 

The availability of the ghost moth host (G) is crucial to the nematode (N); how- 
ever, while the ghost moth's in-year dynamics are modeled explicitly, they are reset 
at the beginning of each wet season, and the model becomes a single-species yearly re- 
turn map for the nematode. This is also the case for the modeled alternate host (H). 
Clearly, the ghost moth has a large and detrimental effect on the lupine, but this 
interaction is not modeled here. The persistence of the nematode is so precarious, 
and the health of the lupine is so dependent on the parasitoid's control of the ghost 
moth, that it suffices to consider the (generally transient) survival of the nematode. 
In Chapter |4| we search for a positive influence on the soil-dwelling nematode's per- 
sistence by examining the possible effects of a second host species. 



Chapter 2 

Spontaneous Patchiness in a 
Host-Parasitoid Integrodifference 
Model 

2.1 Introduction 

Host-parasitoid and other victim-exploiter interactions have long been a staple of 
mathematical investigations in ecology. An interesting subset of such systems are 
those in which the victim's and exploiter's dispersal behavior is crucial to the outcome 
of their interaction. In many cases, the exploiter has significantly greater motility 
than the victim (e.g. [I9]), which can cause the formation of intriguing and even 
counterintuitive spatial patterns [1]. 

When the tussock moth and its parasitoid exploiters, described in Chapter [T| are 
modeled by a pair of one- dimensional reaction-diffusion partial differential equations 
with negligible victim diffusion, a cursory singular perturbation analysis suggests the 
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possibility, and predicts the shape, of a steady state solution [5]. The predicted 
solution consists of an outbreak region of nontrivial host density and a region in 
which only parasitoids exist. Surprisingly, the host density is highest near the edge 
of the outbreak, outside of which it falls rapidly to zero. This result is in qualitative 
agreement with data from the field [22] ■ More recent work has found that such a 
pair of reaction-diffusion equations produces traveling waves with fronts of the shape 
described above - with highest host density near the front |38]. However, that work 
also places fairly strict limits on the types of victim growth behavior and boundary 
conditions that may produce spatially interesting steady state solutions. 

The main drawback to the partial differential equation formulation of the model is 
that it is continuous in time, while host-parasitoid systems often are not. The tussock 
moth is univoltine and it would therefore be more sensible to model the interaction 
in discrete time, as a pair of integrodifference equations. Integrodifference equations 
bring with them the additional benefit of greater flexibility in choosing the kernel 
describing the dispersal behavior |27]. A first step in this direction is taken in |18] . 
where an integrodifference model is built upon the Nicholson-Bailey parasitoid model 
with density-independent host growth. After introducing and analyzing our model 
for the tussock moth system and demonstrating its spatial behavior, we will return 
to that earlier model for comparison. 

We will now formulate our integrodifference model. Following this, we explore the 
dramatic spatial patterns possible for our model, and determine the requirements for 
their existence. First, using a singular perturbation approach we find that an AUee 
effect in the growth of the host is necessary for a patchy final state. We then visualize 
the process leading to the formation of a patchy spatial distribution through numerical 
simulation. With this process in mind, we form a more concrete understanding of 



2.2. The Model 



6 



the behavior of the system en route to its final state, first by relating that behavior 
to the local dynamics of our model, then by linear stability calculations, and finally 
by finding an estimate of the maximum spatial extent of a host outbreak. 



2.2 The Model 



Similar to |1H], the basis for our model will be a Nicholson-Bailey host-parasitoid 
interaction: 

iVi+i = iV,e^(^^)-''^*, (2.2.1) 

Pi+i = Nt{l- e-"^') . (2.2.2) 

The function / provides some form of density-dependent host growth. In this model, 
such density dependence is vital to the possibility of stable coexistence [21]. As 
an example, simple density-dependent (Ricker) growth would be given by f{N) = 
r(l-iV) [2J. 

For any x,y in our problem domain Q, let kd{x,y) be a dispersal kernel - a 
density function for the probability that an organism at y at time t will be at x at 
time t + 1 [29]. The parameter d will describe the magnitude of dispersal, similarly 
to the diffusion coefficient in the continuous-time model. Then our spatial model is 

Nt+^{x)= [ fc,(x,y)iVi(y)e^(^*(^))-'^^'(^)rf2/, (2.2.3) 



Pi+i(x) = / kDix,y)N,iy) (l - e~^''^^^^) dy, (2.2.4) 



n 



where D ~ 1 and e D, both positive, are the dispersal parameters for the parasitoid 
and host, respectively. 



2.3. Singular Perturbation Analysis of Steady State 

In much of the analysis to follow, we will use the Laplace kernel, 



kdix,y) = j^e-\^-y\/''. (2.2.5) 



This dispersal kernel is reasonable for many species p9j and has the additional benefit 
of mathematical tractability, as will be seen. Similar to previous models we will take 
our domain, for most of this paper, to be some interval [0, L] G M. 

2.3 Singular Perturbation Analysis of Steady State 
2.3.1 Regular Perturbation 

To obtain a first-order approximation to a time-invariant solution (A^, P), we formally 
take the limit as £ — * 0. Since kd{x, y) is an approximate identity in the sense of 



convolutions as d 0, (2.2.3) becomes N{x) = A^(a;)e^(^(^))-"-f'(^), so N{x) = 
or e-^*^^*-^^^""'^^^-' = 1. In the latter case, f{N(x)) = aP(x), and we will say that 
N{x) = f^^{aP{x)) in a sense to be explained below. 



From (2.2.4) we see immediately that choosing N(x) = can only lead to the 
trivial solution {N,P) = (0,0). If we choose the other approximation, differentiating 
twice gives 

P" =^{P- r\ciP) (1 - e-^"")) with (2.3.1) 

P'(0) = -ip(O), P'{L) = -^P{L), 

an equivalent problem to the integral equation for P{x) as e ^ 0. Note that P = 
is also a solution to this problem. 



2.3. Singular Perturbation Analysis of Steady State 



2.3.2 Transition Layers 



Following a derivation similar to that of (2.3.1), the integral equation for N{x) with 
£ > is equivalent to the problem 



e^N" = iV (1 - e^(^)-'^^) with 



N'(0) = -N(0), N'(L) = --N{L). 
£ e 



(2.3.2) 



(2.3.3) 



The previous regular approximation would lead us to believe that A^'(O) < and 



N'{L) > 0. This is because, as is evident from (2.3.2)-(2.3.3), there is a layer at each 
of the domain endpoints in which is changing rapidly. These transition layers are 



not surprising given that (2.2.3)-(2.2.4) allows loss at the boundary and the dispersal 



of the host is small. The properties of the solution in and near these layers can 
be determined by singular perturbation analysis, which we omit here. Suffice it to 
report that they are not significantly different from previous results; the host density 
is highest at the edge before declining rapidly near the boundary. 

Since we have found two possible regular approximations, we now turn our atten- 
tion to the existence of internal transition layers. If such layers can exist, it would 
signal the possibility of striking spatial patterns, with stable patchiness - i.e. separate 
coexistence and extinction subdomains - despite the underlying spatial heterogeneity 
of the model. Suppose, then, that a transition layer occurs at some Xq G Q. Because 
the dispersal of P is much higher than that of A^, the density P should be approx- 
imately constant across the transition. So let P ^ Pq there. Define the rescaled 



variable ^ = {x — XQ)/e. Rewriting (2.3.2) in terms of ^, we have 



N" = N{1- e^(^)--^o) 



(2.3.4) 
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Suppose without loss of generality that the transition occurs between a subdomain 
where N{x) = on the left and a subdomain where N{x) = f~^{aP{x)) for P{x) ^ 



on the right. Then the "boundary" conditions for (2.3.4) are 



lim N{0 = and lim N{^) = /"^ (aPo) • 

>— oo ^^+oo 



(2.3.5) 



It is instructive to consider (2.3.4) as a pair of first-order ordinary differential 
equations for and A^'. That system has fixed points at {N,N') = (0,0) and 



(/^^(aPo)j 0). The conditions (2.3.5) require these fixed points to be saddles with a 
heterocline connecting them. In order for them to be saddles, the determinant of the 
Jacobian must be negative at both points: 



det J(0, 0) = e^^o)-"^^" - 1< and 

det J {r' (aPo) ,0)=r' (aPo) /' (/"^ (aPo)) < 0. 

That is, /(O) < aPo and /' (/~^ (c^-Po)) < 0. The first of these implies that, if / is 
nonincreasing on its domain, f{N) = aPo cannot be solved for positive A^. But the 
existence of transition layers requires such a solution. The second inequality implies 
that / is decreasing at some point in its domain. Therefore interior transition layers 
are an impossibility if / is monotonic. 

Suppose that / satisfies the conditions for saddles at the fixed points in question. 
Then there remains the additional necessity of a heterocline connecting them. We 



determine how to satisfy this condition by manipulating (2.3.4) subject to (2.3.5): 



N"N' = N{1- e^(^)~'^^°) N', so 
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\ I {m')'d^ = 1 iV (1 - e^(^)-'^^«) iV'cie, so 

= I N{1- e/W-«^o) ^de, so 



r\aPo) 



j N{1- e^(^)-'^^«) dN = 0. (2.3.6) 



The first fact evident from (2.3.6) is that there must be two values of at which 
f{N) = aPf). If there were only one such value, it would be the only possible definition 
ofN = (aPo), and either f{N) < aPo, or f{N) > aPo, for all G (0, /"^ (aPo))- 
If this were the case, the integral could not vanish. Therefore /^^ (ctPo) must be 
defined as the larger value of A^ at which f{N) = aPo (here we dismiss as biologically 
unlikely a growth function / that crosses aPo at three points). For continuous /, 
since / must be decreasing at /^^ (aPo), it must be increasing at the other solution 
of f{N) = oPq. So / has the unimodal form typical of an Allee effect. 

2.3.3 Summary of Results 



We have analyzed our host-parasitoid model (2.2.3)-(2.2.4), probing the possibility 
of a patchy steady state - a long-term spatial distribution with both extinction and 
coexistence subdomains, sharply segregated. The conclusion of our rigorous analysis 
is that such spatial patterns may only occur if the growth of the host has an Allee 
effect. That is, at small enough host densities, per-capita growth must decrease as 
density decreases. 




Figure 2.3.1: Pattern formation for D = 10, e = 0.1, a = 2, f{N) = {1- N){N -0.2), 
with initial outbreak of width 24 in the left part of the domain. Solid and dashed 
curves represent host and parasitoid densities, respectively. 
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2.4 Numerical Experiments 

Our analytical results thus far have focused exclusively on the final state of the 
system. As will be seen shortly, the dynamics that lead to that state are at least 
as mathematically interesting and biologically important as the long-term spatial 
pattern itself. Moreover, the pattern has been described only in general terms up to 
this point. We now present a visual example of the process of pattern formation in 
the host-parasitoid model. 



2.4.1 Laplace Kernel 



We have explored the behavior of (2.2.3 )-(2.2.4) with the Laplace kernel (2.2.5) 



through extensive numerical simulation, using a fast Fourier transform-based con- 



volution algorithm, as described in p?|, with various grid sizes. Figure 2.3.1 shows 
a typical pattern formation and the steps that lead to it. The initial conditions are 
shown. After oscillating somewhat, by t = 40 the densities inside the outbreak closely 
obey the regular approximation N{x) = f~^{aP{x)). This continues as the outbreak 
spreads until about t = 540, when at the center of the outbreak falls below the 
level predicted by the regular approximation and becomes locally extinct by t = 750. 
Eventually another local extinction occurs and the patches arrange into the pattern 



shown at the bottom of Figure 2.3.1 



Numerical tests with an assortment of parameters and admissible growth functions 
often yield similar results - the initial outbreak spreads until it seems to reach a 
critical width, at which point it divides into two outbreaks, which continue to spread 
and divide until filling the domain. The movement of the outbreak's front is entrained 
by the spread of the host, which is slow and steady, as expected from results in |28] . 
Changing the magnitude of parasitoid dispersal d is roughly equivalent to changing 
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the length of the domain. The dependence of long-term behavior on other factors is 
detailed and analyzed below, and the mathematical arguments are fully in agreement 
with numerical results. 



The calculation of (2.3.6) is analogous to calculations given in ^ and |3H| for sim- 
ilar partial differential equation models. However, analogies to the further conditions 
derived in |38] appear impossible here, given the numerical evidence that interior 
transition layers form in pairs for this model. 



2.4.2 Other Kernels 

Though our rigorous analysis has focused on the Laplace kernel, one of the strengths 
of integrodifference models is their adaptability to the varying dispersal behavior of 
organisms. This variability is reflected in the numerous possible dispersal kernels 



that can be used in (2.2.3 )-(2.2.4), depending on the organisms under study. 



We have simulated the model (2.2.3)-(2.2.4) with a variety of dispersal kernels 



using the parameters given for Figure 2.3.1 The variance of each kernel k(i{x, y) was 



scaled to match the variance of the Laplace kernel (2.2.5) for each d. Figure 2.4.1 



shows the numerical steady state result for the logistic kernel 



kd{x,y) = -^sech^ ( "V^ 
Asd V 2sd 



the Gaussian kernel 



the double gamma kernel 



CTdV^Tr 



-\^-y\/dd 



kd{x,y) = (x - y f 
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Figure 2.4.1: Steady states for other dispersal kernels. 



and the double WeibuU kernel 

For these last two kernels the shape parameter - referred to as a in [31] - is taken to 
be 3, yielding bimodal functions with kd{x,x) = 0. 

The locations of extinction areas seem to depend remarkably little on the specifics 
of the kernel used. Even the double gamma kernel, derived from biological assump- 
tions leading to dramatically different dispersal behavior [3Yj, and also the double 
WeibuU kernel, result in a very similar long-term spatial configuration. The principle 
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difference evident with these and the Gaussian kernel - in short, the non-leptokurtic 
kernels - is that interior patches lack the characteristic increase in host density near 
the edge of outbreaks, noted in 



2.5 Routes to Heterogeneity 



We now attempt to understand why pattern formation in our host-parasitoid model (2.2.3 )- 



(2.2.4) proceeds as described above, through a process of slowly spreading host out- 
breaks repeatedly punctuated by outbreak divisions. To that end we will employ 
arguments of varying mathematical rigor. 

2.5.1 The Nonspatial Model 

Before further consideration of the dynamical behavior of our spatial model, it will 
be useful to bear in mind some of the properties of the underlying nonspatial equa- 



tions (2.2.1 )-(2.2.2). In the N-P plane, there are four nuUclines of the difference 



equations (see Figure 2.5.1). Of particular interest is the curve P = f{N)/a, a 
nullcline for A^, and the curve N = P j [l — e~"^), a nuUcline for P. 

There are three fixed points along the P = axis, and we will make the biologically 
reasonable assumption that there is exactly one more along the nullcline for A^. This 
point lies in the positive quadrant only if a > 1 (we take the carrying capacity of 
to be 1). When a is not much larger than 1, the fixed point is oscillatory but stable. 
As a increases, the nullcline for P is shifted up relative to the nullcline for A^ and the 
fixed point moves to the left along the nullcline for A^ and loses stability in a Hopf 
bifurcation. The other fixed points on the nullcline for A^ are never stable. If the 
AUee effect is strong - i.e., the lesser zero of f{N) is positive - then the fixed point 
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Figure 2.5.1: NuUclines of the nonspatial model. 



at the origin is stable. Otherwise it is not. 



2.5.2 Evolution to the Steady State 

Returning to the spatial model, consider initial conditions given by a single, narrow 
patch of host and parasitoid in the interior of an otherwise empty domain. Inside the 



patch, because of low dispersal the host density approximately obeys (2.2.1). The 



parasitoid, on the other hand, is subject to higher dispersal into the empty part of the 



domain, where it is lost. Referring to Figure 2.5.1 the nuUcline for P is essentially 
shifted downward and the interaction at any point inside the outbreak is stabilized, 
approaching a point on the nullcline for N relatively quickly. 

As the outbreak slowly spreads, the effect of parasitoid dispersal on local dynamics 
in its interior is reduced and the nullcline for P approaches its true position, as 
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determined by the parameter a. The interaction at any point inside the outbreak is 
entrained by the movement of the intersection of the curves, slowly moving left along 
the nuUcline for toward the nonspatial equilibrium. 

If the nonspatial equilibrium lies to the left of the maximum of the nuUcline for 
A^, as the local interaction moves past the maximum it loses the stability temporarily 
imparted by the dispersal of the parasitoid. If the Allee effect is strong, as discussed 
above, the origin is a stable fixed point, and the local dynamics approach it. The 
first point to approach extinction in this way is the center of the outbreak, since the 
parasitoid density is least affected by dispersal loss there. As extinction is approached, 
the parasitoid density is maintained away from zero by an influx from nearby points 
that are not yet approaching extinction. So the host density at the center of the 
outbreak is driven to zero; parasitoids that disperse to the center are lost. In effect, 
two separate outbreaks form, and as they spread the process described above is 
repeated near the center of each of them. This spreading and dividing continues 
until spread is halted at the edges of the domain. 

In Figure |2.5.2 a time series is plotted for the densities at a single point near 



the center of the outbreak shown in Figure |2.3.1[ from t = to 750. 40 time steps 
are required to reach the temporary fixed point on the nuUcline for A^, after which 
the densities move along the curve for 500 time steps until reaching its peak and 
departing for the P axis. 



2.5.3 Dependence on Parameters and Initial Conditions 

As noted above, the formation of a heterogeneous steady state requires a nonspatial 
coexistence fixed point to the left of the maximum in the nullcline for A^ (or equiva- 
lently, to the left of the maximum of f{N)) and a strong Allee effect. With a weak 
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Figure 2.5.2: Time series at outbreak center for Figure 2.3.1 



AUee effect, transient pattern formation is observed during numerical studies, due to 
tfie higli-amplitude oscillatory nature of the map and its slowing near the origin's 
saddle. 

For a reasonably low AUee threshold, a positive coexistence fixed point results 
for large values of a. It is possible to achieve stable spatial patterns for such a, well 
beyond even the range in which there are limit cycles in the nonspatial model. As is 



clear from (2.3.6), a and Pq are inversely proportional for any given /. For Pq? the 



density of the parasitoid at the transitions, to decrease, the fraction of the domain in 
which coexistence occurs must be reduced. This is observed in numerical simulations, 
as coexistence regions become narrower and farther apart with increasing a. In some 
cases, these regions occur far from the domain edges, with complete extinction near 
the boundary. 



Figure 2.5.3 gives a holistic view of the formation of spatial heterogeneity. Pa- 
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Figure 2.5.3: Time series for on the whole domain for parameters and initial 



conditions from Figure 2.3.1 Host density is shown in shades of gray; i.e., black 



represents absence of hosts. 



rameters and initial conditions are as before. Note that the essential features of the 
stable spatial structure are formed within the first 2000 time steps, although the ex- 
tinction subdomains take longer - much longer than we have plotted - to reach their 
final positions. 

The asymptotic behavior of the model may sometimes depend on the initial condi- 
tions. If the initial conditions are extensive enough (for example, a nonzero constant 
throughout the domain), limit cycles can result. Numerical work indicates that, as 
one example, with D = 10, e = 0.1, a = 2.1, and f{N) = (1 - N){N - 0.1), stable 
heterogeneity arises for initial conditions of limited extent while stable limit cycles 
arise for an extensive initial state. 

Certainly for an infinite domain, if the nonspatial dynamics lead to limit cycles, as 
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with the parameters just given, then the spatial resuh for uniform initial conditions 
will be the same limit cycle at each point in the domain. It is not altogether surprising, 
then, that extensive enough contiguous nonzero initial conditions will lead to limit 
cycles with such parameters. This dependence on initial conditions can be understood 



somewhat more rigorously by consideration of the linear stability properties of (2.2.3) 



(2.2.4) 




Wave Number to 



Wave Number w 



Figure 2.5.4: Linear stability where outcome is dependent on initial conditions. 



2.5.4 Linear Stability 



The linear stability analysis of a system of integrodifference equations is explained 
and elaborated in [37], and carried out on a model similar to ours in [IS]. Briefly, the 



Jacobian of (2.2.1 )-(2.2.2) at the coexistence fixed point is 



l + Nf{N) 



-aN 



\ 



where N is the host density at coexistence. The characteristic function of kd{x) 



1 g \x\/d^ ^Yie Laplace convolution kernel, is k, 



2d 



So, much of the behavior 
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of (2.2.3)-(2.2.4) is related to the eigenvalues A of 



KJ 



^ 1 + Nf'{N) -aN ^ 



v 







J 



\ 



1 _ e-/W Nf{N) ^ 



for each wave number uj. Figure 2.5.4 shows the properties of these eigenvalues for 
the parameters given before, for which long-term behavior depends on the extent 



of initial conditions. In the left panel of Figure 2.5.4 , the magnitude of the largest 
eigenvalue (upper solid curve) and the magnitude of its imaginary part (lower left 
solid curve) are plotted. Eigenvalues of K(ci;)J are strictly real for uj greater than 
about 0.1. In the right panel, more detail is given for small uj. In each panel the line 
|A| = 1 is shown for clarity, and in the right panel the cutoff for complex A is shown. 

There are two intervals for uo in which instability is found. Perturbations away 
from the coexistence densities with low cj, or equivalently, long wavelengths, grow 
oscillatorily. Some perturbations with higher a;, or shorter wavelengths, grow mono- 
tonically. We now relate this to the behavior of the model. Extensive initial conditions 
of the kind we have used, viewed as a perturbation from the coexistence densities, 
have a considerable component of long wavelength. Such components grow much 
faster for the parameters under consideration than components of short wavelength. 
But monotonically growing perturbations of short wavelength are precisely the cause 
of extinction subdomains between outbreaks. For these parameters, with extensive 
enough initial conditions, the wild oscillation of the long wavelengths disrupts the 
growth of perturbations with shorter wavelengths predicted by linearizing near the 
coexistence densities. 



Figure 2.5.5 is analogous to Figure 2.5.4 for the parameters that we have used 



in the rest of this paper (see Figure 2.3.1). The maximum magnitude of the real 
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Figure 2.5.5: Linear stability where spatial patterns are robust. 



eigenvalues is comparable to that of the complex ones, allowing spatial patterns to 
form more readily. 

Note that K(co')J — > J as u; — > and thus, not surprisingly, the behavior of 
the model under spatially extended perturbations matches that of the nonspatial 



model (2.2.1 )-(2.2.2). Now, for K(co')J with real eigenvalues, the larger of these, A+, 



is given by 



2A+ = ke{uj)Jii + A;d(cu) J22 + \j {ke{uj)Jii - A;z)(cu) ^22^ + 4A;£(u;)A;d(cu) ^12^21- 

So if J12J21 < 0, which is certainly the case in our victim-exploiter system, we 
have A+ < ks{uj)Jii- Note that ke{uj) 1_ as toe — > 0, so A+ < Ju, and for fixed uj 
this upper bound becomes more accurate as e decreases. Note also that k£,{uj) 0+ 
as uD — s> 00, so the bound becomes more accurate for fixed u as D increases. In 
summary, if 1/D ^ ^ 1/e, then Ju = 1 + Nf'{N) is a tight upper bound on 
real eigenvalues, so since > 0, instabilities leading to pattern formation are only 
possible if f'{N) > 0. This means that, as deduced earlier from simulations, the 
nonspatial coexistence equilibrium must lie to the left of the maximum of f{N). 
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The parameters we have used are such that if ^ 1, then 0.1 = 1/D uj 



1/e = 10. Hence in Figures 2.5.4 and 2.5.5, real eigenvalues are largest near u = 1 



Also, as the nonspatial equilibrium is moved farther left relative to the maximum of 



f{N), as in Figure 2.5.5, these eigenvalues become larger since f{N) is concave. 



2.5.5 A Bound on Patch Radius 



Consider a single small outbreak centered around some point Xc in the interior of the 
domain, far from the boundary, in which local dynamics have reached their temporary 
equilibrium (i.e. they lie on the nuUcline for A^). As the patch spreads it may reach 
a width at which it divides at its center. As discussed above, this occurs when 
the value of A^ at the center maximizes f{N). Let us call this value Amax? and let 

Pma.x /'(A^max)/'^- 



We have from (2.2.4) that 



Pt+i{xc) = j kn{x,, y)Nt{y) (l - e'^^'^^)) dy. 



Turning again to the Laplace kernel, in order for the patch not to divide it must be 
that 



Since Aj = outside the patch, we have 



Xc+R 



Xc—R 



(2.5.1) 



where R is the radius of the patch. Inside the patch, since the kernel is leptokurtic 
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and we are considering the moment at which the patch divides, we will make the 
approximation Nt ^ Nt{xc) = A^max and Pt ^ Pt{xc) = Pmi,^- Then 



1 
2D 



Xc+R 



Xc—R 



x^+R 



^ j e-l^-^l/^iV^ax (1 - e-"^— ) f/y 



Xc—R 



iVma. (1 - e-^^-) (1 - e-«/^) . 



So the patch divides approximately when 



which is when 



R = Dln 



iVmax (1 - e-'^^— ) 



-aPrr, 



Prr 



(2.5.2) 



Note that this approximation was derived for a single non-oscillatory patch in an 
otherwise empty domain. A patch with neighbors should have a somewhat smaller 



radius since the assumption leading to (2.5.1 ) does not apply. This is why it is seen in 



numerical simulations that patches at the boundary are wider than interior patches; 
at the boundary the assumption that A^^ = outside the patch is partially true. The 



approximation (2.5.2) does not hold for patches with sustained internal oscillations 



since, as explained previously, such a patch may be arbitrarily large. 



Figure 2.5.6 demonstrates the accuracy of (2.5.2) by comparing it to the actual 



maximum size obtainable in simulations for various values of a and the Allee thresh- 



old. It turns out that in most cases (2.5.2) should be an upper bound on R, because 
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Approx. 
Actual 





0.15 

Allee Threshold 



Figure 2.5.6: Comparison of approximation and actual maximum radius. On the left, 
f{N) = {1 - N){N - 0.2). On the right, a = 2.4 and / is quadratic and scaled to 
have maximum 0.2. 



Nt (^1 — e "^*) > A^max (l ~ 6 «-Pmax^ when Nt varies more than Pt, as near the out- 



break's center during a division (see Figure 2.3.1 at t = 540). 



2.6 Discussion 



We have observed that our model (2.2.3)-(2.2.4) for a discrete-time host-parasitoid 



system can exhibit stable spatial structure. As in previous and somewhat analogous 
models, the distribution of the host within outbreaks is qualitatively similar to that 
seen in the field, with increased density near the edge [22]. We have shown that this 
pronounced spatial patterning requires a strong Allee effect in the host and fairly 
unstable underlying (nonspatial) dynamics. In such a situation, low host motility 
actually acts as a survival or stabilizing mechanism. The stabilizing influence of 
spatial heterogeneity - as well as the likely local instability of most host-parasitoid 
systems - is discussed in [25] . 

Spatial patterns form through a process of spread and division. Their formation 
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does not greatly depend on the specific dispersal behavior of organisms, other than 
that the dispersal of the host should be comparatively short. In fact, the necessary 



conditions we derived for internal layers in Section 2.3.2 are independent of the para- 
sitoid dispersal kernel. We saw numerically that patterns may form with other forms 
of both kernels, although the distribution within an outbreak does depend on the 
nature of dispersal. This dependence could be grounds for further analytical investi- 



gation. It is also worth noting that the dispersal kernel properties used on page [22 
are not unique to the Laplace kernel. 

Increased parasitoid efficacy results in more sparse patchiness. That is, coexis- 
tence areas become narrower and farther apart as a means of continuing persistence 
of the overall system. For a given set of parameters, a bound on the possible width 
of coexistence patches may be found. Beyond that width, a patch will spontaneously 
divide and the new patches will move apart. In the next chapter we will further 
explore this behavior and its implications. 

While the spatial patterns we have observed and explained are striking and may be 
unfamiliar in the context of spatial host-parasitoid interactions, very similar patterns 
are found for a continuous-time predator-prey model in [.36j- A somewhat similar 
singular perturbation analysis and parameter requirements are given. However, as 
observed in [37j, less tenable biological assumptions are invoked to achieve those 
patterns. Also, the focus of that work is on Turing instability. For the examples 
given in this paper, Turing instability is impossible because when the nonspatial 
model has a stable equilibrium, patterns cannot form, and this is true for most 
reasonable parameters and growth functions. Rather, our examples demonstrate 
diffusion-mediated stability. Lastly, as noted before, our model is far more relevant 
to tussock moths and their enemies than is that of [36J, or other continuous-time 
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predator-prey models such as those in [5] and [38] . 

As mentioned at the outset of this paper, a model similar to ours, but with density- 
independent host growth, is briefly considered in |1H] . The linear stability analysis in 
that work resembles ours, suggesting the possibility of pattern formation. However, 
stable patterns of the kind seen in this paper are impossible there. Firstly, depression 
of the parasitoid nuUcline by dispersal cannot achieve stability of local dynamics in 
that model. Secondly, in that stability analysis there are no real eigenvalues with 
magnitude greater than 1. Numerical simulations of that system show the same 
oscillatory instability that occurs in the underlying nonspatial model. 

Another model related to the tussock moth system is given in [3l]. It is discrete 
not only in time, but also in space, with three patches coupled by dispersal. This 
simple model allows for precise analysis of the conditions required for spatial patterns 
analogous to those described here. The result of that analysis is that the host nullcline 



must must have a "hump," as with an Allee effect, just as we found in Section 2.3.2 



It is fruitful to compare our model's behavior to the dynamics of tussock moths 
in the field. High tussock moth density near the edge of outbreaks [22] suggests 
that dispersal of parasitoids is leptokurtic, with a dispersal kernel shaped roughly 



like that given by (2.2.5). Failure of experimental moth invasions outside preexisting 
outbreaks [32] is consistent with an Allee effect and with numerical simulations of 
our model. Even more solid evidence is given in [20], where it is reported that the 
nonlinear functional response of predatory ants induces an inverse density dependence 
in tussock moths. The details of this are considered in the next chapter. 

Further questions are raised that have not been answered to date by field obser- 
vations. The stability of the nonspatial interaction is of interest. That is, would the 
tussock moth density fall to zero or oscillate if parasitoids were not lost to dispersal? 
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Also, in an ample and unoccupied habitat, how fast do new moth outbreaks spread? 
How sparse are old outbreaks? Is there spontaneous division of large outbreaks in 
the field? 

In addition to the implications of our model discussed above, there is another 



detail that should be carefully noted. As seen in Figure 2.5.3 and discussed elsewhere, 
because of low host motility, outbreaks spread fairly slowly. But once the basic 
spatial pattern is formed, extinction areas travel orders of magnitude more slowly. 
It is unreasonable to expect the final steady state ever to be reached under such 
conditions, especially in the actual tussock moth and lupine system that motivates 
our model, where the habitat is threatened by other organisms and may change 
dramatically in the course of a few years. However, qualitative knowledge of the 
state toward which the system is moving helps to understand and possibly predict 
its transient behavior. 
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Chapter 3 

The Effects of Habitat Quahty on 
Patch Formation 

3.1 Introduction 

The trophic cascade brought about by soil-dwelhng nematodes, whereby the pro- 
ductivity of the lupine is protected to some extent from root-feeding enemies, is an 
important point of interest to the tussock moth system modeled above. We do not 
endeavor to formulate a unified model for the entire food web, from the nematode to 
the tussock moth's enemies (and also including the ants which prey on tussock moth 
pupae), in this work. However, since some of the connections between the various 
components likely only operate in one direction - i.e. the tussock moth has no effect 
on either nematode or ant dynamics - it may suffice to explore the possible back- 
ground conditions, in the form of model parameters and other assumptions. While 
we ignored outside factors in the preceding chapter, certain qualities of the larger 
natural system have a large influence on the striking patterns observed there. 
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We found certain general necessary conditions for pattern formation in Chapter |2] 
We will now see how the feasibility of these conditions is affected by the nature of the 
tussock moth's habitat. We first investigate in detail the shape of the host nuUcline 
given the environmental factors present in the natural system. We then turn to the 
possible results of varying the habitat quality, which determines, in part, the position 
of the host nullcline. 

3.2 An Allee Effect Induced by Saturating Preda- 
tion 

3.2.1 Modeling the Growth of the Host 

Much progress has been made recently in connecting mechanistic reasoning about 
continuous-time in- year dynamics to apparently phenomenological models of discrete- 
time growth. In |l13J, a population consisting of adults and juveniles is considered. 
Depending on the kind of interactions assumed to occur between and within these 
classes in continuous time during the year, various year-to-year maps are derived 
for the density of adults. Likewise, a continuous consumer-resource interaction for 
in-year dynamics is used in \16\ to derive, by varying the specifics of the consumer, 
another assortment of discrete-time models. Similarly, a "semi-discrete" model for 
a host-parasitoid system is considered in [13] ; parasitism is modeled as a continuous 
process throughout the larval stage of the host, and all other processes are modeled 
in discrete time. 

There are a number of natural mechanisms that may lead to growth with an Allee 
effect [3, m]. Some of these mechanisms, such as the difficulty in finding a mate 
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at low densities, have been considered explicitly in discrete-time models |8], |3], Il2] . 
The methodology of using a continuous-time model for in-year dynamics to derive a 
discrete map is extended to Allee growth mechanisms in [12], where various mate- 
finding behaviors are considered, along with cannibalism. 

It has become apparent that predator-prey interactions can induce Allee effects in 
a variety of situations. In one example, a model predator that attacks only a certain 
stage of its prey has been shown to exhibit a growth threshold |16]. As outlined 
in [T3], saturating predation can induce an Allee effect in a prey population. This 
general mechanism is considered in ^2] for a discrete-time growth model. 

In this section, we focus on a population which is victim to stage-specific, general- 
ist, saturating predation. In particular, we model a species with a pupal stage during 
which it is susceptible to such a predator. This describes certain members of the 
family Lymantriidae, such as the tussock moth studied in Chapter [2] and the gypsy 
moth. As mentioned previously, pupae of the tussock moth in coastal California have 
been found to be subject to attack by ants, a generalist predator [20]; pupae of the 
gypsy moth are preyed upon by mice [TT]. We will formulate two models for such a 
victim population, the first an alteration of that given in [42.J, and the second using 
an approach similar to [13], [TB], and [15] . 

We wish to find a map of the form A^j+i = F (Nf), where Nt represents the density 
of the adult female population at the beginning of the winter before year t. These 
adults lay eggs which overwinter and become larvae. The larvae feed and are subject 
to density dependent survival as a consequence of limited resources. The survivors 
pupate; during the pupal stage, they are victim to a generalist, saturating predator. 
The females that survive this stage to become adults comprise the population A^t+i, 
lay eggs, and so forth. 
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The functional response of tlie saturating predator, or tlie rate at whicli it con- 
sumes prey per capita, will be taken as 



(3.2.1) 



1 + sx 

where x is the prey density, m represents the strength of predation (incorporating 
the constant density of predators), and s is related to the handling time or satiation 
of the predator. 

3.2.2 Escape Probability 

A discrete-time model for an AUee effect induced by saturating predation is formu- 



lated in 112]. Based on the functional response (3.2.1 ), the probability of an individual 



in a population with density escaping predation for the entire season is 

m 

/ [Nt) = e i+'^N ^ 

so with Ricker density dependence, the year-to-year map is 

Nt+1 = I {Nt) Nte''^^-^) = Nte^ T^)~^+^t^ (3.2.2) 

where K is the carrying capacity absent predation. This map, however, assumes 
predation throughout the season, simultaneous with density dependence. In the case 
considered in the present paper, resource limitation occurs first: 



Fo = Ar,eni-ifJ, 
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where Yq is the density of larvae that begin the pupal stage. Then the probability of 
escaping predation during pupation is of the form 



e ^+"^0. 



so the year-to-year map is 



3.2.3 Continuous Predation 

The formulation above only evaluates the consumption rate (3.2.1) at the beginning 



of the pupal stage. If, instead, predation is considered as a continuous process during 
the pupal stage, we model it as 



^ mY 



1 + 



where Y is the rate of change of the pupal population. Integrating and applying the 
initial condition, at the end of the pupal stage we have InF + sY = Inlo + sYq — m. 
Here we have either taken the length of the pupal stage to be 1, or equivalently 
rescaled m. This relation between Y and Yq may also be written 

Ye'^ = Foe'^""", (3.2.4) 

in which form we see that the relation defines an increasing function Y (Yq) for positive 



^0- In fact, from (3.2.4) it immediately follows that Y (Iq) = -W (sloe'^°~™), where 



W is the Lambert W function [5j (see also [TH] for a similar ecological application of 
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W). 



Again we model the larval stage with Ricker dynamics: 



The Ricker growth equation has been derived in at least two mechanistic contexts 
- cannibalism [101 [13] and a limited resource [T^. We use it here because it is a 
plausible model for the development of eggs (produced in numbers proportional to 
Nt), through the resource-limited larval stage, to the beginning of pupation. The 
year-to-year map is 



(3.2.5) 



3.2.4 Properties of the Maps 



Recall the form of the Nicholson-Bailey model (2.2.1 )-(2.2.1 ). In Chapter |2| we found 
certain necessary conditions under which the host-parasitoid integrodifference model 



(2.2.3)-(2.2.4) can exhibit dramatic and stable spatial patterns, a typical example of 



which is shown in Figure 2.3.1 To attain these patterns, the host must exhibit a 



strong AUee effect and the coexistence fixed point (see Figure 2.5.1) of the nonspatial 



model (2.2.1 )-(2.2.2) must lie to the left of the maximum in the host's nuUcline. Since 



the host nullcline is P = f{N)/a, the shape of / directly determines the shape of 
the nullcline. For small P, the parasitoid nullcline can be accurately linearized as 
P ^ 2N — 2/a. It crosses the N axis at = 1/a and rises steeply. 

In Chapter |2| as in a model for the gypsy moth in [30], a phenomenological growth 
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Figure 3.2.1: Growth function / for maps (3.2.3) and (3.2.5) witli varying pre- 
dation intensity and liandling time. In each case, the function is plotted for 
K = 0.7, 0.75,. ..,1.0. 
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function of the form 

/ (Nt) = r (1 - Nt/K) {Nt - c) (3.2.6) 

is used to produce a liost population with carrying capacity K and Allee threshold 
c. Conveniently, this function produces a host nuUcline that is symmetric about its 
maximum, providing ample opportunity for the steep parasitoid nuUcline to intersect 
to the left, as required for stable outbreaks. However, this cannot be taken for granted 
in natural systems, given the variety of mechanisms that may shape the left side of 
the nuUcline. 

Many of these mechanisms are reviewed, for example, in [71 01]. Some of the most 
common are those related to the increased individual difficulty of finding a mate at 
very low population levels. A few models of the resulting Allee growth, derived from 
underlying individual mating probability, are reviewed in [3]. These models have 
in common a kind of singular behavior near zero population - exemplified by the 
explicitly boundary-layer model presented and empirically validated in [23] - so that 
/ (Nt) = In {Nf^i/Nt) is extremely steep to the left of its maximum. 

In light of the the results of Chapter |2| a natural means to compare models is 
by the comparison of their growth functions /. Clearly, not all discrete-time single- 
species models are of the form Ni;_^_l = Nte-^^^^\ but they can be rewritten as such. 
For any model Nt+i = G {Nt), we wiU let / (A^*) = In {G {Nt) /Nt). We assume that 
G {Nt) /Nt has a limit as Nt 0, and define / accordingly there. 

Note now that 

dNt+i _ .(0) 



dNt 



Nt=0 



The derivative at Nt = for the map (3.2.3) is ™', just as calculated in |42j for 



the map (3.2.2). We may differentiate the map (3.2.5) at A"^ = using the chain rule 
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and the fact that W'{0) = 1^: 



Nt=0 



We have /(O) = r — m in each case, so each map exhibits a strong Allee effect if 



r < m. Of course the phenomenological map, with / defined in (3.2.6), has a strong 
Allee effect if c> 0. 

It is interesting to note that whether the Allee effect is strong does not depend on 
the carrying capacity K in any of the models. In the phenomenological model, the 
threshold density c is set as a parameter independent of K. In our other models, the 
threshold density should not be expected to depend greatly on K (so dc/dK ^ 0), 
since /(O) does not depend on K. The positive fixed point A^* of each of our models, 



however, moves with K (that is, dN*/dK > 0), as can be seen in Figure 3.2.1 



In our mechanistic models, the strength of the Allee effect depends, in a sense, on 
the predator-related parameters m and s. As we decrease m, predation becomes less 
intense; as we increase s, the predator is more easily saturated and so predation be- 



comes more density-dependent. It is clear that, as m — > 0, the map (3.2.3 ) approaches 



the Ricker map. That is, the function / converges to / (Nf) = r (1 — Nt/K), which is 
Ricker's growth function. This can be seen by comparing Figures 3.1(a) and 3.1(e)[ 
Similarly, as s ^ oo, for Nt ^ 0, f converges to Ricker's growth function. But since 



/(O) = r ~m for the map (3.2.3 ), / does not converge to the value of Ricker's growth 



function at zero. That is, / converges pointwise almost everywhere as s ^ oo, but it 
does not converge uniformly. The slope of / near Nf = grows arbitrarily large as s 
increases. Also, at every other point, / approaches Ricker's growth function - which 
is decreasing in Nt - so the maximum of / moves toward zero as s increases. This is 
seen by comparing Figures 3.1(a) and 3.1(c)[ 
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For the map (3.2.5), since xe^ = W ^{x), for small m we have 



s \ / s ' 



Nt/K) 



so again the map converges to the Ricker map for any A^^. The above also holds for 
large s, if A^^^ 7^ 0. 

The behavior of the models as s increases mirrors the singular behavior of models 
derived from individual mating probability, mentioned above. In the limit of large 
handling time, then, it becomes unlikely that host and parasitoid nuUclines cross in 
the manner shown in Chapter |2] to be necessary for spatial pattern formation. On 
the other hand, this behavior is not encountered as m decreases. So for reasonably 
small s, the growth function / for each of our models is qualitatively similar to the 



simple phenomenological function (3.2.6). The importance of the behavior of the 



mechanistically-derived function / with changes in - roughly the same as for a 
phenomenological, quadratic function - will now be shown. 



3.3 Consequences of Habitat Quality 
3.3.1 Introduction 

The biological system motivating this work, which centers on the bush lupine, con- 
sists of two host-parasitoid interactions, one in the foliage of the lupine (the tussock 
moth subsystem) and one near the roots. As noted earlier, tussock moths have little 
to no effect on the year-to-year survival of the lupine; however, root-feeding ghost 
moth larvae can do great harm to the plant [19]. We will turn our attention to the 
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particulars of the ghost moths and their parasitoids in Chapter |4j 

The consequences of habitat degradation comprise an important problem in the- 
oretical ecology. A change in the quality of a habitat is usually included in models 



such as (2.2.3)-(2.2.4), in which that quality is not dependent on the rest of the sys- 



tem, by varying a parameter representing carrying capacity. One famous example of 
this can be found in [3^. Indeed, the biological system with which we are concerned 
has been so treated in 119], in which the effects of lowering carrying capacity in a 
reaction-diffusion partial differential equation model of the system are investigated. 
It is found in that work that the size of an outbreak is positively related to carrying 
capacity; as the latter increases, outbreaks grow. 



We will now see that the opposite result holds for our model (2.2.3)-(2.2.4) 



3.3.2 Dependence on Carrying Capacity 



As seen in Section 3.2.4, an AUee growth function derived from the mechanism of 
saturating predation shares many features with a simple quadratic function, including 
its behavior upon a change in carrying capacity. Considering this behavior in light of 
Figure 2.5.1 and the discussion and results in Chapter [2| some of the consequences 
of raising or lowering carrying capacity in the model are evident. 



Recall that the P nullcline in Figure 2.5.1 depends only upon the parameter a. 
As noted in the previous section, the nullcline crosses the N axis at = 1/a and 
rises steeply. As K decreases, the maximum of the A^ nullcline moves left and down 



(Figure 3.2.1). If the coexistence point was slightly to the left of that maximum with 



A' = 1 - as it is for Figure |2.3.1| - then with lower K it might be on the right. So 
decreasing K has a stabilizing effect on the underlying, nonspatial dynamics of the 
model, which means that it reduces or even completely removes the patchiness of 
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the spatial model's steady state. Likewise, increasing K destabilizes the underlying 
dynamics and can cause patchiness to appear. 



This result can be seen more exactly in (2.5.2 ), the estimate of maximum outbreak 



radius. For small Pmax, we have 



T-, ' max 
R ~ D— 



A^max -l/a 

where A^max is the location of the maximum of / (and therefore the nullcline). 
If K is lowered such that Nj^.^^ decreases toward approximately 1/a, the estimate 
of outbreak radius increases without bound. In actuality, as A^max decreases toward 



this singular point, the tendency to oscillate (discussed in Section 2.5.4) overcomes 
the tendency to form steady spatial patterns. As K, and consequently A^max, is 
further decreased, the nonspatial dynamics move through the oscillatory parameter 
range and eventually become stable, and likewise the spatial model reaches a smooth 
steady state devoid of patchiness. 

3.3.3 Numerical Examples 



Figure 3.3.1 demonstrates the dependence on carrying capacity explained above. In- 
creasing K leads to more patchiness, with smaller outbreaks, and decreasing K has 
the opposite effect. As Nmax approaches the singular point, the final state becomes 
more sensitive to the value of K, as expected. Note that, for consistency with Chap- 
ter [2} we will continue to use the usual quadratic growth function, since it appears 



from Section |3.2.4| to be qualitatively valid. 

Another interesting question about the effects of habitat quality is whether the 
overall abundance of the host, measured by N{x)dx (where N{x) is the steady 
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Figure 3.3.1: Consequences of varying carrying capacity, with D — 10, e — 0.1, 
a = 2, f(N) = (1 - N/K){N - 0.2), with initial outbreak in the leftmost quarter 
of the domain. Solid and dashed curves represent host and parasitoid densities, 
respectively. Each simulation was run io t — 10000. 
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Figure 3.3.2: Consequences of varying carrying capacity on overall host abundance 
in a domain of length 80. 



state), increases with K. It would be natural to expect this, were it not for the 



dependence of patchiness on K described above. Figure 3.3.2 shows the numerical 



results on the domain Vt = [0,80], with parameters from Figure 3.3.1 and K varying 
from 0.5 to 1.1. Results are obtained by running the simulation for a given value 
of K until ||A^t(a;) — Nt^iQQ{x)\\i2 falls below some threshold, then slightly varying 
K and repeating the process on the previously-obtained steady state. The evident 
trend is that abundance increases with except near K = 0.9 where two extinction 
regions form. Simulations have not been carried out in enough detail to comment on 
the exact nature of changes in abundance during the transition to patchiness, and as 
noted above and in the next section, a steady state may not always be possible in 
this parameter range. However, it is clear that, while increasing habitat quality will 
lead to greater host abundance for most values of K, this is not always true. 
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3.3.4 The Paradox of Enrichment 

The consequences of varying K, the carrying capacity of the host in the absence of 
the parasitoid, are somewhat counterintuitive but should be familiar to mathematical 
ecologists. They are a reflection of the so-called paradox of enrichment [H] whereby 
a victim-exploiter interaction is destabilized by enriching the habitat of the victim - 
or, in a model, by increasing its carrying capacity. 



This "paradox" has been studied for nonspatial models such as (2.2.1 )-(2. 2.2 ). 
The addition of spatial dispersal, however, has a profound effect on the stability of 
the system under changes in carrying capacity: 



spatial 



/ stable 


lim. eye. 


patchy \ 


^ stable 


limit cycles 


extinct 



nonspatial 



The horizontal line represents carrying capacity K. As it increases, the nonspatial 
model loses stability, first to limit cycles and then to extinction. The extent of the 
oscillatory parameter range is considerable; this diagram was generated for K from 0.7 
to 1.1. The spatial portion of the diagram was produced with the parameters used in 



Figure |3.3.1[ on a domain of length 160, with simulations run to 10000 time steps and 
categorized by eye - the values of K at boundaries between differing behaviors were 
not found with any precision and therefore are not labeled. Rather, the importance 
of this diagrammatic overview is qualitative. Spatial considerations slightly extend 
the parameter range that produces a long-term stable (not patchy) solution. The 
limit cycle range is greatly reduced in size, and beyond it the formation of patches 
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prevents global extinction, which numerical simulations show to be the case up to at 
least K = 2. 



3.4 Consequences of Habitat Heterogeneities 
3.4.1 Introduction 

A more general situation than the previous discussion of habitat quality is the pos- 
sibility that carrying capacity K depends on location x. As we will see, the most 
interesting habitat heterogeneities are abrupt changes, or discontinuities, in carrying 
capacity. 

Discontinuities may be a reasonable expectation given the nature of the habitat 
being modeled. As will be discussed in the next chapter, though the bush lupine 
habitat may be spatially continuous from the perspective of the tussock moth, it 
consists of individual plants with taproots set at some distance from each other. The 
health of each plant depends on conditions in its rhizosphere such as the dynamics 
of detrimental ghost moth larvae and the parasitic nematodes that exploit them. It 
is evident that if there is any dynamical coupling between rhizospheres, it is at best 
weak and sporadic [TU] ; as such, even neighboring bushes may differ greatly in quality. 
We now consider the consequences of such heterogeneities. 



3.4.2 Patch Formation 



As in Section 2.5.5 consider a point Xc inside the domain of the model (2.2.3) 



(2.2.4) near which the local dynamics are settled to the outer approximation N(x) 



f ^{aP{x)) (see Section 2.3.1). Recall that local extinction occurs when the para- 



sitoid density reaches the maximum of the host nuUcline, due to - in a homogeneous 
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environment - an outbreak growing large enough to sustain that density of parasitoids 
inside. 

In a heterogeneous environment, if the point under consideration is in a region 
of low quality (low carrying capacity) relative to nearby habitat, there may be in- 
teresting repercussions for patch formation. Lowered carrying capacity results in a 
lower host nuUcline, which is more easily overcome by the parasitoid density. In a ho- 
mogeneous habitat this is countered by the movement of the coexistence fixed point 
of the local dynamics, which for low enough carrying capacity makes it impossible 
for parasitoid levels to reach the maximum of the host nuUcline. However, a nearby 
region of higher quality habitat can provide levels of parasitoid influx sufficient, when 
combined with locally-generated parasitoid densities, to cause local extinction. While 
homogeneously low-quality habitats are less likely to become patchy through the for- 
mation of local extinctions, in a heterogeneous environment low-quality regions are 
the most likely locations for extinctions. 

More interesting effects can be seen if habitat heterogeneities are fairly abrupt. 
For example, if a domain is divided into one high- and one low-quality region, and the 
drop in carrying capacity is sufficiently discrete and severe, an outbreak beginning 
in the high-quality habitat may be halted at the heterogeneity. For a completely 
discrete (stepwise) drop in carrying capacity, there is clearly a sufficient condition on 
the magnitude of the drop required to halt the outbreak. At the edge of the outbreak, 
parasitoid density has some finite value because D, the dispersal parameter for the 
parasitoid, is relatively large. If the low carrying capacity is such that the maximum 



of the resulting host nuUcline (Pmax in Section 2.5.5) is less than that density, the 
outbreak should not be able to spread past the location of the drop in habitat quality. 
In fact, numerical analyses, which we will now present, suggest that the discrete drop 
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need not be even that severe. 



3.4.3 Numerical Examples 




Figure 3.4.1: Simulation with i^' = 1 in the left half of the domain and K = 0.9 in 



the right half. Parameters and initial outbreak size are the same as in Figure 2.3.1 



We investigate the consequences of a habitat heterogeneity by dividing the domain 



from Section 2.4 into two (equal) regions of higher and lower K. The former will have 



K = 1. Setting K = 0.9 in the low-quality region, the numerical results are shown 
in Figure [3.4.1 An outbreak beginning in the high-i^ subdomain fills it in the usual 
way and spreads into the area of lower K. A local extinction occurs at the sharp 
boundary between the regions, on the side of lower K. The approximate final state 
is shown. 

Similarly, an outbreak beginning in the lower-quality region spreads to the in- 
terface and continues into the subdomain with K = 1. Shortly thereafter, local 
extinction occurs at the interface as above, and the final state is identical. 

Further lowering the carrying capacity in the low-quality region leads to the effect 
described in the previous section - an outbreak beginning in the high-capacity area 
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Figure 3.4.2: Simulation with = 1 in the left half of the domain and K = 0.8 in 



the right half. Parameters and initial outbreak size are the same as in Figure 2.3.1 
In the bottom right panel, the initial outbreak was placed on the left. 



may not be able to spread beyond it, even though spread in the opposite direction is 



still very much possible. Figure 3.4.2 shows spread in the direction of increasing K 
when K = 0.8 in the lower-quality region. It is qualitatively similar to the case with 
K = 0.9 there. Also shown is the apparently final state - the result is identical for 
widely varying grid resolutions - resulting when the initial outbreak is placed in the 
higher-quality subdomain. Note that the outbreak does not divide, because it cannot 
grow quite wide enough. In the previous cases, division was made possible by influx 
from the low-quality region. 

As in Section |3.3 we explore the effects of varying the habitat quality - now 
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Figure 3.4.3: Consequences of varying K2 on overall host abundance in a domain of 
length 80. 



heterogeneous - on the overall abundance of the host. As above, we divide the 
domain into two (equal) regions, with K = Ki in one and K = K2 in the other. 



In Figure 3.4.3, we set Ki = 0.75 constant and vary K2 from 0.5 to 1.1 (the same 



range as in Figure 3.3.2). As K2 approaches Ki, the overall abundance of the host 



increases; however, as K2 grows larger, the abundance begins to decrease, until K2 ~ 
1. Throughout this parameter range, the steady state has a single extinction region, 
on the low side of the heterogeneity as usual. But near K2 = I, a second division 
occurs in the high-i^ subdomain, so that the steady state is similar to that seen in 



Figure |3.4.2[ Contrary to the result for homogeneous habitat quality, the overall 
abundance of the host actually increases when this division occurs. 

The two counterintuitive behaviors evident here - a negative relationship between 
carrying capacity and abundance, and an increase in overall abundance when a local 
extinction occurs - are directly related. An abundance of hosts in the high-quality 
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subdomain drives the formation of an extensive extinction region in the other subdo- 
main by exporting parasitoids. When the productive subdomain experiences a second 
local extinction, that portion of the subdomain is unproductive with respect to the 
highly dispersive parasitoids, and the first extinction region narrows accordingly. 



3.5 Extension to Two Dimensions 

We will now briefly see that all of the qualitative results in one spatial dimension, 
explored in Chapters [2] and |3} apply analogously in two dimensions. We use the 



model ( 2.2.3 )-( 2.2.4) with Laplace kernel 



kdix,y) = ^^e~\\^~y\y''. (3.5.1) 



We simulate the model with a two-dimensional fast Fourier transform convolution 
algorithm (very similar to that in Chapter [2]) on a square domain. In all of our 
simulations we use f{N) = (1 — N/K){1 — 0.2), a = 2, D = 10, e = 0.3, and domain 
n = [0,240] X [0,240]. 

In our first three examples we compare the outcome of a small, random initial 



outbreak on a domain with homogeneous K of various values. Figure 3.5.1 shows the 
initial condition used in all three cases, along with some snapshots of the spread of 
the host with K = 1. As the outbreak spreads, a local extinction occurs in its center. 
Later, a ring-shaped extinction begins to form. However, the initial irregularities in 
the outbreak's edged lead to inclusions which disrupt its symmetry. Eventually a 
second near-ring forms, and spread continues to the boundary. With a circular initial 
condition, extinctions would occur in concentric rings until the outbreak reached the 
domain edge and disturbances propagated back into the domain, breaking symmetry. 




Figure 3.5.1: Two-dimensional simulation with K = 1: Initial condition, t = 450, 
t = 680, t = 1250, t = 2500. Host density shown in shades of grey from = 
(black) to = 1 (white). 
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Figure 3.5.2: Two-dimensional simulation with K = 0.8 (left) and K = 1.2 (right) 
ai t = 2500. Host density shown in shades of grey from = (black) to > 1 
(white) . 



Figure 3.5.2 shows the effects of raising and lowering carrying capacity on the 
whole domain. With K = 0.8, the outbreak fills the domain with no extinction 
areas. With K = 1.2, the long-term state is patchier. 

To simulate heterogeneity, we place a roughly circular region with i^' = 1 in the 
center of the domain, and lower K outside this area. In Figure 3.5.3[ K = 0.9 in 
the lower-quality region. Note that this value of K is high enough for patchiness in 
a homogeneous environment. We place an initial outbreak in the upper-right of the 
domain, completely outside the region with higher K. The outbreak begins to form 
extinction areas, and easily spreads into the higher-quality habitat. Local extinctions 
begin to form around the edges of that area, and then inside it. Finally, the domain is 
filled and local extinctions largely group around the edge of the region with higher K, 
though unmistakably on its exterior (the region is visible as a rough circle of greater 
host density). 

If the initial outbreak lies inside the region of higher carrying capacity, it easily 
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Figure 3.5.3: Two-dimensional simulation with heterogeneous domain quality {K = 1 
in the center, K = 0.9 at the edges): t = 450, t = 750, t = 1500, t = 4000. Host 
density shown in shades of grey from N = (black) to = 1 (white). 



spreads outside and the long-term state is qualitatively identical to that in Fig- 
3.5.3[ If we lower K outside the high-quality region to K = 0.8, which as we saw 



ure 



above produces no patchiness if the habitat is homogeneous, the outcome is essen- 
tially the same. Hosts can spread across the heterogeneity in either direction, and in 
the end there is a solid ring of local extinction around the high-quality region, and 
now with no local extinctions away from it. However, if we further lower carrying 
capacity outside to i^T = 0.6, an outbreak beginning inside the region of high carrying 
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Figure 3.5.4: Two-dimensional simulation with heterogeneous domain quality {K = 1 
in the center, K = 0.6 at the edges): t = 2000, t = 4000. Host density shown in 
shades of grey from = (black) to = 1 (white). 




Figure 3.5.5: Two-dimensional simulation with heterogeneous domain quality {K = 1 
in the center, K = 0.6 at the edges) at t = 6000. Host density shown in shades of 
grey from = (black) to = 1 (white). 
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Figure 3.5.6: Two-dimensional simulation with K = 1 and a = 4: t = 
500, 1000, 1500,2000,4000,6000. Host density shown in shades of grey from = 
(black) to = 1 (white). 
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capacity cannot escape it. This is shown in Figure [3.5. 4[ The high-capacity region is 



soon filled, but spread is stalled at its edge and the interior extinctions simply begin 



to rearrange. On the other hand, as shown in Figure 3.5.5 , if the outbreak starts 
from the low-quality region, it can spread throughout the domain. 

Returning briefiy to the homogeneous case with K = 1 throughout the domain, 
if we take a to be quite large, the shape of the host's spread becomes interesting 
and similar to that observed for a partial differential equation model with an AUee 



effect in the victim. Figure 3.5.6 shows the resulting simulation with the same initial 



condition as in Figure [3. 5.1 The host's spatial behavior - coalescing into thin ribbons 



and breaking into patches while spreading - is reminiscent of the spread of the prey 
in [Sn] , where it is termed "patchy invasion" . 



3.6 Discussion 

We began this chapter with an analysis of an AUee effect induced by saturating pre- 
dation. In many victim-exploiter models, the exploiter nuUcline is a vertical line. 
One prominent example, for a continuous-time model, is found in [H]. In such cases, 
the specifics of the host nullcline - other than its humped shape - may not be im- 



portant. However, with the nonspatial model (2.2.1)-(2.2.2), in which the parasitoid 
nullcline has finite slope, our implicit assumptions require that the host nullcline have 
relatively small slope to the left of its maximum, as it does with the growth function 
we used in Chapter [2] 

While the quadratic function we have used to describe growth is not realistic for 
many instances of the Allee effect in nature, we saw that it might closely match 
the properties of a growth function derived from mechanistic principles based on the 
reality of the system motivating our model. Specifically, an Allee effect induced by 
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saturating predation is likely to generate a growth function of very similar shape and 
behavior to our quadratic approximation. 

We saw numerically how increased predation on the pupal stage of the host 
strengthens the AUee effect. Importantly, and in contrast to Allee growth not in- 
duced by a predator, the host nullcline is roughly quadratic even for a relatively 



weak Allee effect (Figure 3.2.1 for low m). We found the parameter condition under 
which the effect is strong - i.e., such that the host is bistable in the absence of par- 
asitoids. Evidence from the field suggests that the effect is indeed strong there pU] . 
Finally, we saw how the induced growth function, and therefore the host nullcline, 
changes with a change in carrying capacity K. 

The movement of the growth function upon a change in K has interesting impli- 
cations for the effects of varying habitat quality. This is largely due to the paradox of 
enrichment [H], but its destabilizing effects are mitigated by dispersal, with spatial 
pattern formation striking a balance between stability and extinction. Also, in the 
cases we considered numerically, increasing carrying capacity within any parameter 
range that produces a fixed number of patches seems to lead to increased overall host 
abundance, though the abundance falls slightly upon patch division. 

Although spatial spread is not stopped in a homogeneous environment for our 
continuous-space model, as it may be for metapopulation models with homogeneous 
sites [171 [26], the heterogeneities required to stop spread are fairly mild. We have also 
seen how adjacent regions with differing habitat quality affect one another. Similar 
situations, for continuous-time models, are outlined in [H]. The effects of varying 
the quality of only a portion of the habitat can be counterintuitive. If the carrying 
capacity in one area is made sufficiently large, the overall abundance of the host may 
actually decrease. This result is conceptually similar to theoretical studies of marine 
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systems which have predicted that, in many cases, estabhshment of a reserve could 
lower the overall abundance of the protected species in the presence of predation [35] 
or infectious disease [HH] . 
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Chapter 4 

Persistence of Parasitic Nematodes 
Augmented by a Scarce Alternate 
Host 

4.1 Introduction 

As explained in Chapter [T| the roots of the bush lupine provide shelter and sus- 
tenance to ghost moth larvae, which are in turn parasitized by entomopathogenic 
nematodes. In this way, the nematodes act to promote the health of the shrubs in a 
trophic cascade. However, a single infected ghost moth larva can produce hundreds 
of thousands of nematodes, leading to pronounced year-to-year population cycles, 
which can result in the local extinction of the nematode when stochasticity at low 
numbers is taken into account [9]. 

The ability of nematodes to disperse between roots and form a metapopulation is 
questionable, which calls into question the global persistence seen in nematode-ghost 
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moth systems in the field [10] • The aim of this chapter is to explore a new possibility 
to explain the observed persistence of nematodes - the infiuence of an alternate host. 

Cases in which the presence of a second host is pivotal to the survival of a par- 
asitic population, even if that host is inferior to the first, have been observed both 
theoretically [24J and in the field To model a second host for the nematode, we 
will build upon the deterministic model from [9], in which wet season dynamics (from 
t = to t = T) are given by 

H'{t) = -kuHit) - (3H{t)N{t), H{0) = Hq (4.1.1) 

N'{t) = -kNN{t) - (3H{t)N{t) + (3A{t - T)H{t - T)N{t - r), (4.1.2) 
with the next year's initial nematode numbers, A^(0), determined by 

T 

XoN{T) + \iP j A{t)H{t)N{t)dt. (4.1.3) 

T-T 

First we will examine the dynamics of this original model in slightly greater detail 
than previously done. We will then formulate a new model accounting for a second 
host. Finally, through numerical simulation and reasoning from the original model's 
dynamics, we will see how the addition of another host, even in very small numbers, 
may positively infiuence the persistence of the parasitoid. 

4.2 Details from the Original Model 

We will see shortly how certain properties of the model given in [9] may contribute to 
the ability of an alternate host to enhance the persistence of the parasitoids. It will 



4.2. Details from the Original Model 



60 



be fruitful to examine these before moving on to the formulation of the new model. 
4.2.1 In- Year Dynamics 

The behavior of the model during an overexploited wet season is vitally important 
for our purposes since this behavior leads to the dangerously low densities of the 
subsequent season. In truth, the danger caused by overexploitation of the host is not 
the removal of the host from the system, since it is replenished every wet season. The 
danger lies in the premature removal of the host, such that no infections occur late 
enough in the wet season for the resulting nematodes to be protected from the high 
mortality of the dry season. 




Days Days 



Figure 4.2.1: Wet season dynamics during the overexploited year of the two cycle 
with /3 = IQ-^ (left) and /3 = 2 x 10"^ (right). The threshold time T - r = 125 is 
marked with a vertical line. Other parameters: = 0.0001, /ctv = 0.063, Ao = 10~^, 
Xi = 10-3, A(t) = min{10,000e° °9*, 800,000}, Hq = 100. 

Recall that infections occurring after the time T — t produce nematodes with 
reduced mortality during the dry season [9J . For parameter values producing violent 
two-cycles in the yearly map, the last infections in the overexploited year take place 
before this crucial window. However, it is important to note that these infections 
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may occur late enough that the nematode density during the window is quite high. 



Figure |4.2.1| demonstrates these dynamics for two values of the infectivity (3. For 
(3 = 10^"^, the nematodes reach their greatest numbers at the beginning of the window. 
For lower f3 this peak occurs even later. Even with (3 = 2 x 10"'^, nematode numbers 
are significantly higher in the crucial window than in much of the rest of the season. 

4.2.2 Bifurcation Diagram of the Map 

Let us take note of the location and nature of the fixed point and two-cycles of 



the yearly map. Figure 4.2.2 gives a bifurcation diagram on the parameter range 
producing the stable two-cycle that motivates this work. For much of this range, the 
fixed point is at a fairly small density. Notice also that as the infectivity (3 increases, 
the fixed point becomes stable through a flip bifurcation before the two-cycle vanishes 
in a fold bifurcation. The basin of attraction of the fixed point becomes quite large, 
and extends downward almost to = 1. 



4.3 The New Model 

4.3.1 Unsuitability of Deterministic Modeling 

During the wet season, we consider the infection of the alternate host i^ait with a 
per-host infection rate proportional to parasitoid density A^. Each infection produces 
Aait nematodes after some incubation period, which we will take to be the same as in 
the primary host, r. Over the course of the dry season, all soil-dwelling nematodes 
experience high mortality Aq. Nematodes remaining in alternate hosts at the begin- 
ning of the dry season, like those in primary hosts, experience lower mortality, which 
we also take to be A,. 
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Figure 4.2.2: Bifurcation diagram for values of P producing a two-cycle. Solid curves 
indicate stable fixed points or cycles, dashed curves unstable. Parameters are as given 



in Figure 4.2.1 



Written in the deterministic (mean field) form of (4.1.1 )-(4. 1.3), our model would 



be 



H'{t) = -kHH{t) - (3H{t)N{t), H{0) = Ho 
Kitit) = -(3.itN{t)H,,,{t), H,,t{0) = H, 



N'{t) = -kNN{t) - (3H{t)N{t) + (3A{t - T)H{t - T)N{t - r) 



-P^ltN{t)H,lt{t) + /?altAalti^alt(t - T)N{t - t). 



with the next year's initial nematode numbers determined by 



T T 

XoN{T) + K\f3 J Ait)H{t)N{t)dt + /S^itA^it J H^,,{t)N{t)dt 

T-T T-T 
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However, if an alternate host for the nematode exists, it is likely not nearly as 
abundant as the primary host. As such, we wish to examine the repercussions of very 
few alternate host individuals per year on nematode persistence, and so modeling 
this secondary host in a deterministic fashion is out of the question. 



4.3.2 The Model and Its Simulation 

The alternate host and its interaction with the nematode population will be modeled 
as a continuous-time stochastic process, with each individual host subject to an infec- 
tion rate of f^^itN. In the event of an infection, three transitions occur: the alternate 
host population is reduced by 1, the nematode population is reduced by 1, and r 
days after the infection, Aait new nematodes are produced - unless the dry season 
has begun. 

The stochastic process described above occurs simultaneously with the determin- 
istic wet season dynamics. This is numerically simulated by calculating, at each 
time step of the differential equation solver, the approximate probability that a given 
alternate host will be infected during that time step: Ps,itN(t)At, where At is the 
length of the time step. Events are then carried out according to simulation of the 
corresponding random variable for each alternate host. The next year's initial nema- 



tode density is given by (4.1.3), with the addition of the term AjAait multiplied by 
the number of infections that occurred after the cutoff time T — r. 

In all simulations we will use a fourth-order Runge-Kutta routine with step size 



At = 0.1, and the parameters from Figure 4.2.1, unless otherwise indicated. Limited 



tests with varying At indicate that this method is stable and accurate. 
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With the model formulated, we will now see how dynamics such as those described in 



Section 4^ contribute to enhanced persistence in the presence of a scarce alternate 
host. 

4.4.1 Mitigated Crashing 

The obvious means by which the addition of an alternate host to the highly cyclic 
victim-exploiter interaction can enhance nematode persistence is by increasing nema- 
tode densities in low years. 

The number of infective juveniles present in the soil at the end of an overexploited 
season, in the crucial window during which an infection may produce new nematodes 
sheltered from the high mortality of the dry season, can be orders of magnitude 



greater than most of the rest of the season, as seen in Figure 4.2.1 So, given a 
proportional relationship between nematode density and infection probability, rare 
alternate host infections may be most likely to occur during this window. 



Figure |4.4.1| demonstrates the dependence of persistence on /?ait, as well as the 
results of increasing the number of alternate hosts available in each wet season. We 
have used (3 = 2 x 10"'' and Aait = 25, 000, averaging on 20,000 randomized trials 
at 100 points along the /?ait axis. In each simulation we start with = 5 and 
count persistence until falls below 5. Mean persistence time seems to increase 
proportionally with the number of alternate hosts available each year. Also note, 
most strikingly, that there is effectively a /5ait threshold above which the alternate 
host has no effect on persistence. 
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Figure 4.4.1: Dependence of persistence on alternate host infectivity and number of 
individual alternate hosts per year. 



4.4.2 Transient Survival 

An alternate host in our cyclic interaction can enhance nematode persistence by 
raising nematode density to an intermediate level near the fixed point of the year- 
to-year dynamics, perhaps producing a lengthy coexistence transient. As discussed 
previously, the fixed point, even when unstable, often occurs at a relatively low 



density. For the parameters used in Figure 4.2.1 with /3 = 2 x 10~^, the fixed point 
is around N = 200. So if an infection late in the overexploited wet season were to 
produce around 200,000 sheltered infective juveniles, by our model parameters there 
would remain just slightly more than 200 nematodes at the beginning of the next wet 



season. The resulting cobweb diagram is shown in Figure 4.4.2 In this particular 
case, nematodes do not return to low numbers (less than 10) until five years after the 
late-season infection. 
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Figure 4.4.2: Cobweb diagram of transient dynamics after a fortuitous infection at 
the end of an overexploited wet season. 



Figure 4.4.3 shows the results of increasing the productivity of alternate host in- 
fections. Note that doubling Aait to 50,000 has little to no effect on persistence, but as 
AaitAj approaches the fixed point, persistence is increased. However, infection output 
must be increased eightfold to have roughly the same effect as tripling the number 
of individual alternate hosts. This disparity is due to the fact that only a single 
alternate host infection in the crucial window is necessary to lengthen persistence by 
two years. The probability of such an infection increases roughly proportionally with 
additional hosts; increasing productivity has no bearing on that probability, and the 
transient effects of a given alternate infection are only lengthened significantly by 
choosing unreasonably specific nematode densities. 
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Figure 4.4.3: Dependence of persistence on alternate host productivity. 



4.4.3 Deterministic Stability 



As seen in Figure 4.2.2 for larger values of the infectivity /3, the yearly map exhibits 
not only a fixed point at low density, but stability of that fixed point with a fairly large 
basin of attraction. For such values, a single alternate infection before an otherwise 
low-density year could move the dynamics into asymptotically stable coexistence. 



Figure 4.4.4 shows some of the dramatic results possible when infectivity of the 
primary host is chosen such that the system has a stable fixed point in addition to its 
limit cycle. In this case, (3 = 4.0 x 10"'', a value close to the flip bifurcation (thus the 
fixed point's basin of attraction is relatively small). We have also set /3ait = 10~^° and 
Aait = 25, 000, with one alternate host individual available per year. The simulation 
is started on its limit cycle. 

Meaningful alternate host infections during overexploited years of the two-cycle 
may be highly unlikely, since peak in-season nematode densities occur well before 
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Years 



Figure 4.4.4: Time series of dynamics with two attractors. 



the crucial window for such a high value of (3 (see Figure 4.2.1); however, once an 
infection occurs, dynamics stay near the fixed point - and away from nematode 
extinction - for extended periods. Note that the dynamics approach the fixed point 
and are perturbed away by alternate infections frequently. Meaningful infections take 
place much more often in this range of densities than in the overexploited years of 
the limit cycle, because here nematode density peaks later in the wet season. In 
fact, these infections can take place often enough to move the system out of the 
basin of attraction, allowing it to approach the dangerous limit cycle. For example, 



this happens at 2,210 and at 3,730 years in Figure 4.4.4, Note, however, that at 
3,200 years the system moves out of the basin of attraction, only to be immediately 
returned to it by an alternate host infection. 



Clearly, the length of time represented by Figure 4.4.4 is not meant to be realistic 
for a time series of the dynamics in a natural rhizosphere; it simply serves to demon- 
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strate the possible persistence of each of the stable states of the original model under 
the influence of an alternate host. 

4.5 Discussion 

The efficacy of our alternate host in lengthening the persistence of the nematode is 
very dependent on the value of f3, the infectivity of the primary host. However, the 
means by which persistence may be extended are such that, for much of the range 
of values of (3 that produce violent two-cycles, the presence of the alternate host can 
have dramatic effects. For low (3, meaningful alternate host infections - those which 
produce new nematodes sheltered from the high mortality of the dry season - become 
quite likely during overexploited years. For high /3, meaningful infections may become 
less likely, but only a single infection is needed to approach stable persistence. 

It is crucial to the survival of the parasitoid population that some host organisms 
be infected at the end of the wet season. We have seen how the nature of the wet 
season dynamics allows for this if a host exists that is significantly less susceptible 
to the nematode. This host, even in small numbers, then maintains the parasitoid 
during a crucial period when the primary host is unavailable. This is reminiscent 
of an alternate prey of predatory mites, studied in [17], which is present earlier in 
the season than the prey which these mites are meant to control, thus allowing the 
predator population to grow larger - and therefore more effective - before the targeted 
species' numbers rise significantly. 

Note that our conclusions would hold even if the only difference between the pri- 
mary and alternate hosts were their infectivities, (3 and /5ait; respectively. That is, 
the model and results would apply to a single host population in which a few indi- 
viduals are significantly less susceptible to parasitism, through genetic disposition. 
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an effective refuge, or some other means. In such a situation we can see that, para- 
doxically, the existence of individuals that are much less susceptible than the rest of 
the population should have a very positive effect on nematode persistence. 
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